DCA"*" : Dynamical Cluster Approximation with continuous lattice self-energy. 
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The dynamical cluster approximation (DCA) is a systematic extension beyond the single site 
approximation in dynamical mean field theory (DMFT), to include spatially non-local correlations 
in quantum many-body simulations of strongly correlated systems. We extend the DCA with a 
continuous lattice self-energy in oder to achieve better convergence with cluster size. The new 
method, which we call DCA"*" , cures the cluster shape dependence problems of the DCA, without 
suffering from causality violations of previous attempts to interpolate the cluster self-energy. A 
practical approach based on standard inference techniques is given to deduce the continuous lattice 
self-energy from an interpolated cluster self-energy. We study the pseudogap region of a hole-doped 
two-dimensional Hubbard model and find that in the DCA^ algorithm, the self-energy and pseudo- 
gap temperature T* converge monotonously with cluster size. Introduction of a continuous lattice 
self-energy eliminates artificial long-rage correlations and thus significantly reduces the sign problem 
of the quantum Monte Carlo cluster solver in the DCA+ algorithm compared to the normal DCA. 
Simulations with much larger cluster sizes thus become feasible, which, along with the improved 
convergence in cluster size, raises hope that precise extrapolations to the exact infinite cluster size 
limit can be reached for other physical quantities as well. 



INTRODUCTION: 



The study of interacting electrons in a crystalline solid 
remains one of the most challenging problems of con- 
densed matter physics. On a purely theoretical level, 
these models give us insight on spontaneous symmetry 
breaking, which leads to new ground states with exciting 
properties such as superconductivity. On a more practi- 
cal level, the lattice-models allow us to better understand 
materials in which the correlations between electrons de- 
termines its physical properties. The most famous exam- 
ples of such materials are the high- Tc cuprates^ and the 
recently discovered pnictidea^l. A better understanding of 
how the Cooper pairs are formed in these materials might 
lead us in the future to the creation of new materials with 
higher superconducting transition temperature. 

One of the methods of choice to investigate interact- 
ing electrons on a lattice-model is the DMFT'^, which, 
in conjunction with model parameters derived from first 
principles electronic structure calculationPnS, is now ca- 
pable of predicting spectral properties of transition metal 
oxidetP and heavy fermion material^i^Hls). The study of 
the paring mechanism in superconductors, however, re- 
quires inclusions of dynamic correlations between lattice 
sites, and hence the extension of DMFT beyond the sin- 
gle site approximation. To this end, several quantum 
cluster extensions to DMFT have been developed during 
the past fifteen yearpHH. Among these is the DCA^iaHIlI 
a systematic extension to DMFT that includes non-local 
correlations through coarse-graining in momentum space. 
The DCA relies on the assumption that the self-energy 
function is a localized function in real space. In infinite 
dimensions, it has been proven that the self-energy S is a 
delta-function in real space^Sl, in which case this assump- 
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FIG. 1. The positions of the cluster momenta {K} and shape 
of the patches for two 16 site DCA-clusters. Notice that the 
16B site cluster does not have the same point-group symmetry 
as the Brillouin-zone, leading to a lattice self-energy with a 
lower symmetry. 



tion trivially holds. In practice, we see that the locality 
increases with increasing dimension. This explains why 
a single-site DMFT approach generally works better for 
3D materials, but fails to describe materials of quasi ID 
or 2D nature. 

The Dynamical Cluster Approximation was developed 
to study materials of 2D nature, by allowing the self- 
energy to be non-local. In the DCA, the infinite lattice- 
problem is reduced to a finite size quantum cluster impu- 
rity with periodic boundary conditions, embedded into 
a self-consistent mean-field. This reduction is achieved 
via a coarse-graining procedure of the Green's function, 
in which the Brillouin zone is divided into Nc patches 
and the self-energy S is assumed to be constant on these 
patches. In this way, all correlations within the cluster 



are dealt with exactly, while long-range correlations out- 
side the cluster are described via a mean-field. If the 
cluster impurity problem is solved exactly, such as with 
Quantum Monte Carlo (QMC) integration, the DCA will 
reproduce the exact solution of the lattice model in the 
limit of infinite cluster size. 

In practice, the fermionic sign piohleul--'-- imposes an 
upper-bound to the cluster-size and a lower bound to the 
temperature which can be accessed. While small clusters 
have proven to give us an excellent qualitative insight 
on the physical phenomena^^, most physical quantities, 
such as the superconducting transition temperature Tc, 
converge poorly on the available small clusters^S. The 
DCA can therefore not be used as a reliable method for 
quantitive predictions of those observables. 

There are two important factors that influence the 
results of the DCA, both related to the choice of the 
cluster. The most obvious factor is the mean field ap- 
proximation, which reduces the momentum anisotropy 
of the self-energy as the clusters become smaller. One 
can only avoid this error by considering clusters with a 
sufficiently large size. In practice, the critical cluster-size 
is obtained by comparing physical quantities on different 
cluster-sizes. More complicated is the influence of the ge- 
ometry of the cluster. There is a set of different clusters, 
all of which have the same cluster size but different shape 
and therefore different positions of the cluster momentum 
points. In Fig. [l] we show two 16 site clusters for which 
this is the case. The different positioning of the cluster 
momentum points in these two clusters leads to a differ- 
ent geometric shape of the coarse-graining patches and 
thus a different parametrization of the self-energy. This is 
illustrated in Fig. [2] where the momentum dependence of 
the DCA self-energy at the lowest Matsubara frequency 
is shown for the 16A and 16B site cluster introduced in 
Fig. [l] The relative error between the self-energies on 
the different clusters is close to 100% around the Fermi- 
surface, making it unsuitable to derive any quantitative 
results from this calculation. 

One can argue that the influence of the mean-field ap- 
proximation for clusters with the same size is similar. 
Therefore, the difference in results can be brought back to 
the shape of the coarse-graining patches. One example is 
the difference in superconducting transition temperature 
Tc between the 16A and 16B site cluster-^. The role of 
the geometry has been studied intensively by investigat- 
ing the evolution of the magnetic and supercondu cting 
transition temperatures over different cluster sizes^^^^ 
or by comparing the site-oc cupan cies of different clusters 
over a wide range of dopin^^IMi 

The geometric shape dependence of the self-energy is 
built into the DCA by construction, since the DCA self- 
energy is expanded on the coarse-grain patches as^S 

I](fc, tI7„) = ^ 0^^ (fc) Sjj^ {zu„r). (1) 

i 

Here, the set of patches {(pg. (k)} is formally defined 



through the cluster-momenta {Ki}, 



<f>Kik) 



1 Vj : \k- K,\ < \k-Kj 
3j : \k-Ki\ > \k-Kj 



(2) 



and E^ (wm) is the cluster self-energy for momentum 

K^■ 

In this paper, we present an extension to the DCA 
that allows the self-energy to be expanded in an arbi- 
trary large set of smooth basis-functions, and thereby 
itself becoming a smooth function of momentum. The 
inclusion of a smooth self-energy into the framework of 
the DCA requires a new fundamental look at the algo- 
rithm. The resulting extended algorithm will be called 
DCA"*" , indicating an incremental generalization to the 
well-known DCA algorithm. The distinguishing feature 
of the DCA+ algorithm that sets it apart from the DCA 
algorithm is that cluster and lattice self-energies are in 
general different. In the DCA, the lattice self-energy T,{k) 
is a simple extension of the cluster self-energy S^^ via 
the step function form in Eq. M. It therefore has jump 
discontinuities between the patches. In the DCA+ the 
lattice self-energy is a function with continuous momen- 
tum dependence, which, when coarse grained is equal to 
the cluster self-energy. 

The focus of this paper is threefold. First, 

we will present the theoretical background of the 
DCA+ algorithm, without going into any practical de- 
tails. Next, we introduce a practical implementation for 
the DCA+ algorithm and discuss in detail the numerical 
aspects of the lattice mapping implementation. Finally, 
we apply the DCA+ algorithm to the single band Hub- 
bard model in order to investigate the pseudogap behav- 
ior, which has rec ently been investigated in a systematic 
way with the DCA'23211. in the theory section, we first de- 
rive the coarse-graining equations for the DCA+ , which 
define how the lattice system is mapped onto an effec- 
tive cluster problem. This will introduce the key con- 
cepts of the DCA+ approach on a general level. Next, 
we discuss the structure of the DCA+ algorithm in more 
detail. Here, we will pay special attention to the lattice- 
mapping, i.e. the inversion of the coarse-graining, where 
a lattice self-energy is estimated from a given cluster 
self-energy. We will show that the lattice mapping is 
only possible if the DCA assumption of a localized self- 
energy in real space is uphold. In the implementation 
part, we will discuss the lattice-mapping in detail on a 
practical level. In this paper, we propose to perform the 
lattice mapping in two steps. First, we interpolate the 
self-energy obtained from the cluster solver. Next, we 
deconvolute the interpolated cluster self-energy, where 
the patch 0g(fc) is used as the convolution kernel. In 
the physics section, we will use the DCA+ to investigate 
the pseudogap behavior in the low-doping region of the 
two-dimensional Hubbard model. The self-energy in this 
phase is known to be strongly momentum dependent and 
we will show that the pseudogap transition temperature 
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FIG. 2. Momentum-dependence of the DCA and DCA^ self-energies (red and black represent the imaginary and real part) 
calculated on the 16^ and 16B clusters in a half- filled Hubbard model with nearest neighbor hopping f = 1, Coulomb interaction 
U/t = 7 and next-nearest neighbor hopping t' /t — —0.15 at a temperature T — 0.2. For the DCA, one clearly sees a large 
difference between the self-energies of the two clusters at the section (7r,0) — >■ (0,7r), which is close to the Fermi-surface and 
thus physically the most relevant part of the self-energy. In the DCA""" , the self-energies of the two clusters agree very well. 



T* converges faster with regard to the cluster size in the 
DCA+ than in the DCA. 

All calculations in this paper were performed for a 
single-band Hubbard model 



H ^J^^^J^aCja 



UE 



ni^n.ii . 



(3) 



Furthermore, all single-particle functions defined on the 
impurity-cluster are represented by a subscript on the 
cluster-momenta (e.g the cluster self-energy S^), while 
the continuous lattice single-particle functions will have 
the usual dependence on the momentum vector k (e.g. 
the lattice self-energy E(fc) ). An overline over the quan- 
tity signifies that the latter has been coarsegrained. 



J 



Here c]^ (cio-) creates (destroys) an electron with spin 



(7 on lattice site i and 



'^ia'^ia 



is the correspond- 



ing number operator. The hopping matrix tij includes 
nearest (t — 1) and next-nearest (t') neighbor hopping 
and U is the on-site Coulomb repulsion. The effective 
cluster problem of the DCA and DCA+ is solved with 

an continuous-time auxiliary-field quantum Monte Carlo 
algorithnPlSll. 



A. DCA and DCA+ formalisms: 

A system of interacting electrons on a lattice is gener- 
ally described by a Hamiltonian H = Hq + TJjnt, where 
the kinetic energy Hq is quadratic in the fermion opera- 
tors and the interaction Hi^t is quartic. It's free energy 
n may be written in terms of the exact single-particle 
Green's function G as 



I. THEORY 

In this section, we present the generic structure of 
the DCA+ algorithm, without going into any implemen- 
tation details. First, we introduce the key features of 
the DCA+ algorithm that distinguish it from the DCA, 
and show that the latter is just a specialization of the 
former. Next, we present a geometric interpretation of 
the DCA+ algorithm in terms of the functional repre- 
sentation space of the self-energy. This interpretation 
provides guidance for how cluster-dependent features are 
incorporated into the lattice self-energy, and offers in- 
sights for the derivation of a practical implementation 
of the DCA"*" algorithm that will be discussed in the fol- 
lowing section. In order to keep the notation simple, we 
will omit the frequency parameter zu in all equations. 



n[G] = Trln(-G) + $[G] - Tr[(Go"^ - G-^)G] . (4) 

Here we have used a matrix notation for the Green's func- 
tion G of the interacting system described by H and the 
Green's function Go of the non-interacting system de- 
scribed by Hq. $[G] is the Luttinger-Ward functionaP2l 
given by the sum of all vacuum to vacuum "skeleton" dia- 
grams drawn with G. The self-energy E is obtaine d from 
the functional derivative of $[G] with respect to (J^SM 



6^G] 
SG 



(5) 



and is related to the Green's function via the Dyson equa- 
tion 



Gn 



G- 



(6) 



These two relations imply that the free energy is sta- 
tionary with respect to G, i.e. 5Vt\G]/5G = 0. In prin- 
ciple, the exact Green's function G and self-energy E 
can be determined from the self-consistent solution of 
Eqs. ([5]) and ^. However, since the functional $[G] 
is usually unknown, an approximation is required that 
replaces the exact $[G] by a known or a computable 
functional. Conserving approximations replace the ex- 
act ^[G] by an approximate functional, which sums up 
certain subclasses of diagrams that are thought to cap- 
ture the dominant physics. In general, this results in a 
weak coupling approximation. A different approach is 
taken in the DCA: rather than approximating the Lut- 
tinger Ward <&, the functional representation space of 
the Green's function is reduced by replacing the exact 
Green's function G(fc) by a coarse-grained Green's func- 
tion G^ in momentum space defined as 



G^ 



dk(f>g{k)Gik). 



(7) 



where the coarse-graining functions 4>j^{k) have been de- 
fined in Eq. (pi). We note that approximating G in this way 
corresponds to an approximation of the Laue function, 
^fc +fc fc +fc ' which expresses momentum conservation 
at each vertex in the diagrams defining ^MH]^ Yov the 
single site DMFT approximation {Nc = 1), (t>{k) is con- 
stant over the entire Brillouin zone, and consequently the 
Lauc function is replaced by Admft = 1, i-e. momentum 
conservation is disregarded. For a finite size DGA cluster 
{Nc > 1), the Laue function restores momentum conser- 
vation for the cluster momenta K and reads in terms of 
the (j)j^{k) 

ADCA(fcl,fc2,fc3,fc4) = %,+;?3,;?2+i?^ (8) 

X <Pk, (^l) 'I^K2 (^2) (t>K:> (^3) 0i?,(fc4) ■ 

By replacing the exact Laue function with its DCA ap- 
proximation in the Luttinger Ward functional, the mo- 
mentum integrals over the Green's functions in the di- 
agrams defining the ^-functional are reduced to sums 
over the finite set of coarse-grained Green's functions de- 
fined in Eq. ([7|. This way, $[G] becomes identical to the 
Luttinger- Ward functional of a finite size cluster and the 
computation of the corresponding self-energy 



E^-p^ = S<!>[G^]/SG^ 



K 



K 



(9) 



becomes feasible. As such, within the DCA approxima- 
tion the free energy functional Q\G] becomes 



VL 



DCA 



[G] = Trln(-G) + $[G] - Tr[(Go ^ - G-^)G] 



From stationarity of the free energy, 5Q\G]/5G 
one obtains the Dyson equation within the DCA 



G^\k)-G-\k) = Y^ 



^KCk)^T^ 



0, 



(10) 



Here, the right hand side follows from 6Gj^/SG ~ 4>i^{k) 
and 5^Gji\/5Gg- = TPS^^. Eqs. ([tI, (|9l and ([l0| form 
a closed set of equations which is solved iterativciy un- 
til self-consistency is reached. This is the DCA algo- 
rithm. Following Eq. (10), the self-energy S(fc) of the 
lattice Green's function G(fc), which is used to compute 
the coarse-grained Green's function in Eq. ([7|, is approxi- 
mated by a piecewise constant continuation of the cluster 
self-energy S^^ , which changes between different mo- 
mentum patches but is constant within a given patch. 



E(fc) = 5]s^cA^^(fc). 



(11) 



K 



With the DCA+ algorithm we introduce in this paper, 
the DCA framework is extended to allow for a more gen- 
eral relationship between the lattice self-energy E(fc) and 



cluster self-energy Sjj than that in Eq. (11). In the 



DCA"*" , in analogy with Eq. h\, we only"3emand the 
cluster self-energy to be equal to the coarse-grained lat- 
tice self-energy. 



^K 



dk4)g{k)^{k) 



(12) 



In the DCA algorithm, thi s re quirement is trivially sat- 
isfied since according to Eq. ( 11 ), E(fc) is set to the cluster 
self-energy E(-ftr) for momenta k in patch Pi. However, 
it is important to realize that Eq. ( [l2| allows for a more 
general approximation of the lattice E(fc), which, for ex- 
ample, can retain its smooth momentum dependence in- 
stead of the DCA step function character. To proceed, it 
is convenient for our purposes to express the free energy 
as a functi onal of the self-energy. By following the work 
of PotthofpSES^ .^g eliminate the Green's function G in 
favor of the self-energy E to write the free energy as a 
functional of the self-energy S, 



nm = -TTln[-{Go' - E)] + (/:$)[E] . (13) 

Here, the functional (£$)[S] is obtained from <I>[G] 
through a Legendre-transformation 



(/:$)[S] =$-Tr[SG]. 



(14) 



Replacing S(fc) in (£<&)[!]] with the coarse-grained self- 
energy in Eq. (12), i.e. T,{k) sa X^i? <^a'(^)^k' then yields 



(£$)[!]] =$ 



2^ ^ii^K' 



(15) 



K 



where G^ is the coarse-grained Green's function defined 
in E g. ^ . If this functional is used in the free energy in 
Eq. (|13[), one obtains at stationarity, (5ri[I]]/(5S = 0, 



[Go-i(fc)-S(fcri=5]^^(^)G^ 



(16) 



K 



Here, the right hand side foUows from Sflj^/SI] = 

(figik) and (£$)[i;^]/(5i;^ = ~^k- Using the identity 

J dk(j)j^{k)(j)j^,{k) — Sg- ^, and multiplying both sides 

with J dk<j)j^{k) results in the DC A+ coarse-graining 
equation 



Here, aj are the expansion coefficients of the lattice self- 
energy corresponding to the basis- function Bj{k). Con- 
trary to the DCA, the coarse-graining patches (f>f}{k) in 
the DCA+ are not linked in any shape or form to the basis 
functions in which we expand the lattice self-energy. As 
was mentioned in the previous section, the DCA+ maps 
the full lattice problem into a cluster impurity prob- 
lem embedded into a mean field by coarse-graining both 
the lattice self-energy and lattice Green's function. The 
cluster-mapping in the DCA+ is thus very similar to the 
cluster-mapping in the DCA, with the exception that we 
use a continuous lattice self-energy in the coarse-graining 
of the Green's function 



Vbz Jbz 



dk(j)j^{k)T.{k), 



(19) 



G 



Nc 



K 



Vb 



z Jbz 



dk<j)j^{k) [G"{k)]'^ ~Y.{k) 



Gg- 



dkc^ji{k)[G^\k)-nk)\ 



Eq. ( 19 ) can now be simplified by using the explicit 



(17) expansion of the lattice self-energy in Eq. (181 



We note that in contrast to the DCA algorithm, 
the lattice self-energy S(fc) enters in the coarse-graining 
step. It is related to the cluster self-energy S^^ through 
Eq. ( |l2| ), i.e. its coarse-grained result must be equal to 
'^{K). The special choice S(A:) = ^^0^(fc)I]^ satis- 
fies this requirement and recovers the DCA al gori thm. 
But in general, S(A:) needs to only satisfy Eq. (12 1, i.e. 



one has more freedom in determining a lattice self-energy 
E(/c) from the cluster Y,{K). In the DCA+ algorithm, we 
take advantage of this freedom to derive a S(A:) that re- 
tains a smooth fc-dependence and thus is more physical 
than the piecewise constant S(fc) of the DCA. As in the 
DCA, the cluster self-energy S^ may be determined from 
the solution of an effective cluster problem described by 
(£$)[!]] as a functional of the coarse-grained propagator 
Y.[K] = i:[G{K)]. This, together with Eqs. ^ and ^ 
form the basis of the DCA"*" algorithm. 

A detailed description of the algorithm will be given in 
the implementation section. Evidently, determining the 
lattice self-energy E(fc) from the cluster self-energy E^^ 
through inversion or deconvolution of Eq. (12 1 presents 
a difficult task. 



B. Structure of a DCA cluster-calculation: 

Since the lattice self-energy S(fc) no longer is restricted 
to Eq. (fTl), it can be expanded into an arbitrary set of 
smooth basis functions {,Bi(fc)}, such as cubic splines or 
crystal harmonics, i.e. 



E(fc) = ^e,(ft)a,, 



(18) 



%=E 



dk(l)j^^{k)Bj{k) 



(20) 



= Pi 



Here, Pij is a projection operator, defined by coarse- 
graining the basis function Bj over patch i. Note that in 
the DCA, this projection operator becomes the identity- 
operation (5i,j. Hence, the coarse graining of the lattice 
self-energy in the DCA is an implicit operation (ci = 
fi^ ), while in the DCA+ it becomes explicit. 

With the in trod uction of the cluster-mapping in the 
DCA+ in Eq. (20), the lattice mapping is conceptually 



well defined as long as the inverse of the projection- 
operator P exists. Assuming that P~^ exists, we can re- 
trieve the expansion coefficients of the lattice self-energy 
from the self-energy of the cluster-solver S^ in a straight- 
forward manner 



=j:(p-')^ 



^R. 



(21) 



This closes the DCA+ iteration and allows us to carry 
out a self-consistent calculation. 

In Fig.lSJ we have summarized the generic structure of 
the DC A^ algorithm, without specifying yet any imple- 
mentation details of the lattice- mapping. In the " cluster- 
mapping" step, the lattice Green's function and self- 
energy are coarse-grained onto the patches defined by 
$^(A:) to give Gg and E^, respectively. A cluster solver 
algorithm such as QMC is then used to calculate, from 
the corresponding bare Green's function G^ ^, the inter- 
acting Green's function and self-energy E^ on the clus- 
ter. In the "lattice-mapping step", which is missing in 
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menta K. A smooth interpolation will almost certainly 
fail to satisfy Eq. (12), i.e. the main requirement of the 



DCA+ that the coarse-grained lattice S(fc) is equal to the 
cluster E^. Such a procedure was proven in Ref.^ to 
lead to causality violations when the cluster self-energy 
is added back to the inverse coarse-grained propagator in 
the " cluster exclusion" step to avoid overcounting of self- 
energy diagrams. In the DCA'^ , the lattice self-energy 
is different from an interpolated cluster self-energy and 
the self-energy that enters the cluster exclusion step is 
given by the coarse-grained lattice self-energy. Because 
of this, the proof given in Ref.'^Sl does not apply and the 
DCA+ algorithm is not automatically plagued by causal- 
ity problems. Although we do not have a rigorous proof 
that the DCA"*" algorithm remains causal, we have never 
encountered any causality violations in the application of 
this method to the single-band Hubbard model. 



C. Analysis of the projection operator Pi,j and its 
connection to the locahty of E(fc): 



FIG. 3. The generic structure of a self-consistent DCA 
algorithm, in which the cluster- and lattice-mapping play a 
central role in order to connect the continuous lattice self- 
energy E(fc) with the cluster self-energy E^j. Convergence is 
reached when the cluster-solver produces a cluster self-energy 
E^ equal to the coarse-grained self-energy E^ = E (-/?). 



the standard DCA algorithm, a new estimate for the lat- 
tice self-energy S(A;) is then computed through inversion 
of the projection operator Pij. The lattice self-energy 
then enters the next cluster-mapping step via the lattice 
Green's function G{k). In the implementation section 
of this paper, we will describe in detail how the lattice- 
mapping can be done in a numerically stable way. 

Due to the distinction between the lattice and clus- 
ter self-energy in the DCA+ algorithm, we can not use 
the convergence criteria of the DCA. In the latter, con- 
vergence is reached if the self-energy (lattice or clus- 
ter) of the previous iteration is equal to the current 
one. If one monitors only convergence on the lattice self- 
energy in the DCA+ algorithm, one might stop the iter- 
ations although the cluster solver still produces a clus- 
ter self-energy S^ that differs from the coarse-grained 
lattice self-energy f^{K). This would indicate that the 
DCA+ does not converge to a stationary point of the 
free energy functional fi. To avoid such a problem, 
we demand that convergence is reached only when the 
coarse-grained lattice self-energy T,{K) and the cluster 
self-energy 'Sg agree to within the Monte Carlo sampling 
error. 

It is important to note that the proposed algorithm 
is fundamentally different from a simple interpolation 
of the cluster self-energy E^^ between the cluster mo- 



From the previous subsection it follows that the pro- 
jection operator Pij plays a central role in the imple- 
mentation of the DCA"*" algorithm. In order to obtain 
a self-consistent algorithm, it is conceptually clear that 
the projection operator has to be invertible. In prac- 
tice, however, this might not be so easily achieved. In 
this section, we give the reader an intuitive understand- 
ing for this operator and show that its inverse exists if 
the DCA locality assumption is satisfied for the lattice 
self-energy. Furthermore, we discuss how the projection 
operator Pij is influenced by the choice of the cluster. 

To this end, we expand the lattice self-energy in terms 
of cubic Hermite splines^. These functions form a basis 
for cubic splines and obey a convolution property. The 
lattice self-energy can therefore be written as sum over a 
very fine mesh {fc^} in momentum space. 



E(fc) = Y^ cr^^ nik - ki) with E(fc,) 



^fe. 



(22) 



It has to be stressed that choosing Hermite splines as 
a basis will not influence the conclusions we obtain here 
and thus does not reduce the generality of our arguments. 
It just simplifies the discussion, since the expansion index 
i can now be identified with a lattice momentum ki in the 
fine lattice mesh and the expansion coefficient Ui with the 
lattice self-energy at that lattice momentum ki. Next, we 
generalize the cluster-mapping in Eq. ( 20 1 , by replacing 



the cluster momentum points {Ki} by the fine lattice 
{ki}. The coarse-graining then becomes a convolution of 
the lattice self-energy with the patches and we obtain 



1.0 r 



^ 



hli._. 



0.8 



0.6 



0.4 



0.2 



0.0, 



» * 






. . N,, = 4 

. . N,, = 8 

. . N, = 16 (A) 

o o AT,, = 16 {B) 

. . iV,, = 32 

o o N,. = 64 

. . iV^ = 100 



^ *^" »mi rt n r . »^ vl f^ 



H 
D 

hi 

c 
CO 
DC 



50 



40 



A = 16{A), B = 4: 
A = IG{A), B = 8 
A = 16{A), B = 1C){B) 



10 



15 



20 



25 



30 




FIG. 4. The leading eigenvalues of various clusters on a fine 
mesh of 512 points. We can clearly observe a strong decay 
of the leading eigenvalues for small clusters, which becomes 
weaker with increasing the cluster-size. This observation ex- 
plains the intuitive notion that large clusters can describe 
finer features in the self-energy, since the image-space of larger 
clusters contains more eigenvectors. 



FIG. 5. The dimension of the union image space X^. UX^- 
for two different clusters A and B versus the eigenvalue index 
i. Since the rank of I;^. and X^. both equal i, any deviation 
of the rank for the space X^ . U X^. from i indicates that the 
projection operators of clusters A and B span different image 
spaces. One can clearly observe that the differentiation of the 
16A site cluster eigenspace with smaller clusters occurs faster. 



ki / J kj 



dk Mk - h) nik - kj) . (23) 



The projection-matrix Pj j has now become a sym- 
metric, square matrix. The latter allows us to do a 
spectral decomposition of Pj: j: into its eigenspace. If 
we represent its eigenvalues by A and its corresponding 
eigenvector by e\, we obtain 



% = Y. % 51 ^ ^^(^») ^ ^A (fcj) 



(24) 



In terms of the eigenspace of the projection-operator, 
the cluster- and lattice-mapping can now be written as 

cluster-mapping: Sj^ ^ X! "^ ^'^fc ' '^a(^j)) ex{ki) 

A 

lattice-mapping: ctj:^ = ^ A~^ (E^^, eA(fcj)) ex{h) 

(25) 

Here, the inner-product (a, b) is represented by a sim- 
ple dot-product between the two vectors a and b. From 
Eqs. (25 1, it is clear that the spectrum {A} of the 
projection-operator Py plays a central role in the cluster- 
and lattice mapping. In Fig.[4j we show the leading eigen- 
values (i.e. having the largest absolute value) of Pi,j for 



various clusters. One can clearly observe that all eigen- 
values are smaller or equal than one and decay rapidly 
for small clusters [Nc < 8) and slowly for large clusters 
{Nc > 32). This can be easily understood from the form- 
factor of the patches. The latter are very similar to box- 
car filters, which are one of the most common low-pass 
filters used in the field of signal processing. Since the 



coarse-graining of the lattice self-energy in Eq. ( 23 1 can 



be rewritten as a convolution with the patches, the pro- 
jection operator Pij will in fact reduce all the Fourier 
components during the convolution, insuring that the 
i2-norm of any function in the eigenspace never grows. 
Consequently, this is also true for all eigenvectors, which 
leads us to conclude that the eigenvalues have to be less 
or equal to 1. 

Using the spectral decomposition of the projection- 
matrix, we can split the representation space of the con- 
tinuous lattice self-energy into the image-space I and the 
kernel-space /C of the projection operator Ptj. Since our 
projection-operator does not follow the strict mathemati- 
cal definition of a projection operator^, we define the im- 
age Xg as the space spanned by the eigenvectors that have 
an eigenvalue larger than e. Here, e is a small, positive 
cut-off parameter. The kernel /Cg contains the remainder 
of the space, and is thus spanned by the eigenvectors with 
an eigenvalue smaller than e. Due to the inversion of the 
eigenvalue in Eq. (25), the lattice-mapping is only well- 



defined on the image-space Xg. This brings us to the first 
important observation. In order to do a self-consistent 
DCA+ calculation, the coarse-grained lattice self-energy 
should always be entirely defined on the image-space X^ 
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FIG. 6. The correlation between the magnitude of the lead- 
ing eigenvalue and the delocalization of the its corresponding 
eigenvector for various clusters. 



of our projection operator. Otherwise, there exists no 
well-defined transformation that maps the cluster self- 
energy back into the lattice self-energy, which in turn 
breaks the DCA+ self-consistency loop. Notice that this 
requirement holds trivially in the case of the traditional 
DCA, since in that case the projection matrix is simply 
the identity-matrix of size Nc, and all eigenvalues are 
equal to one. 

Eq. (24 1 can also explain how the geometry of the 
patches will influence the results obtained with the 
DCA+ . In Fig. [sl we plot the union space of the image 
spaces I^. and I^, versus eigenvalue index i for differ- 
ent clusters. The plot shows very clearly that the first 



leading eigenvectors are equal to each other, and gradu- 
ally diverge as eigenvectors with smaller eigenvalues are 
added. This brings us to the second observation. If one 
wants to carry out a DCA-calculation with results that 
are independent of cluster shape, the cluster self-energy 
has to be representable on the intersection of the image- 
spaces Xg of both clusters. 

So far, we have only discussed and introduced strict 
geometrical criteria on the lattice and cluster self-energy, 
that indicate when a DCA+ cluster calculation is feasible. 
In order to link geometrical criteria to physics, we show 
in Fig. [6] the delocalization of the leading eigenvectors 
(r^). Formally, we define the delocalization as 



{r')x 






(26) 



At close inspection, we can see a clear correlation be- 
tween the absolute value of the leading eigenvalues A and 
the delocalization of its corresponding eigenvector for all 
cluster sizes. This correlation shows that the space X^ 



is actually spanned by the eigenvectors with a small de- 
localization. As a result, satisfying the geometric crite- 
ria to do a self-consistent DCA+ calculation is essentially 
equivalent to satisfying the DCA-assumption of locality 
for the lattice self-energy. Another important conclusion 
that can be drawn from Fig. |6] is that the number of 
vectors that span the space Xg^o.25 becomes larger with 
increasing cluster size. This correlation reflects the intu- 
itive notion in the DCA that larger clusters can describe 
finer features of the lattice self-energy. 



D. Role of the cluster in the DCA+ 

In the DCA algorithm, the real space cluster takes a 
central role. It completely defines the basis-functions in 
which the self-energy is expanded. Furthermore, the real 
space cluster dictates how the lattice is mapped on the 
cluster through the coarse-graining procedure. Conse- 
quently, solutions obtained with the DCA algorithm usu- 
ally dependent on the particular choice (shape) of the 
cluster. In practice, this leads to a very good qualita- 
tive description of the physics, but prohibits quantita- 
tive analysis, as calculated physical quantities strongly 
depend on cluster shape. In the DCA"*" , we start from an 
expansion of the self-energy into an arbitrary set of basis- 
functions. In this way, the influence of the real space clus- 
ter is reduced, since it does not dictate the basis-functions 
on which the self-energy is expanded. The real space 
cluster only specifies how the cluster is mapped on the 
lattice through the shape of the coarse-graining patches. 
Consequently, the focus in the DCA+ shifts from the real 
space cluster to the projection operator Pij- This oper- 
ator embodies the quantum cluster approximation of the 
DCA"*" , since it connects the cluster self-energy with the 
lattice self-energy in a purely geometric way. The projec- 
tion operator is only defined by the set of basis-functions 
of the lattice self-energy and the real space cluster and 
not subjected in any way to physical parameters (such 
as temperature, band-structure, interaction terms, ...). 
This purely 'geometric' property of the projection oper- 
ator allows us to find a priori the necessary conditions 
to which the cluster self-energy has to be subjected, in 
order to allow for a self-consistent, cluster- independent 
DCA+ calculation. These necessary conditions that fol- 
low from the discussion in teh previous subsections are: 

• In order to perform a self-consistent 
DCA+ calculation, the cluster self-energy has 
to converge in the image-space I^ of the projector. 

• In order to perform a cluster-independent 
DCA+ calculation on the cluster A and B, the 
cluster self-energy needs to converge on the in- 
tersection of the image-spaces of both projectors 
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FIG. 7. The decay of E^ for & N^^ 100-site cluster with 
I] jt — 7 and t'/f = for various temperatures at half filling. 
For high temperatures {T < 0.3), the system is only weakly 
correlated and E_r will rapidly decay. For low temperatures, 
the correlations exceed the cluster-radius Re = 5. 



cation in the localized nature of Wannier orbitals, from 
which the tight-binding models are derived. Since the 
self-energy is a correction to the band-structure due the 
interaction between the electrons, the Wannier interpo- 
lation method seems a suitable interpolation algorithm. 
Okamata et al.^^ have examined this possibility implic- 
itly, by expanding the lattice self-energy S(fc) into the 
cubic-harmonic basis- functions {C^(/c)}. 



I](fc) = ^C^(fc)Ex 



(27) 



K 



CKCk)-^j:-^' 



{K~k) 



This approach only works when the self-energy S^^ is 
sufhciently smooth, such that the real-space self-energy 
YiB converges on the cluster in real space. Notice that 



the latter is implicitly computed in Eq. (27), since 



II. IMPLEMENTATION 

In the last section, we have introduced a projection 
operator P^j and shown its involvement in the cluster 
and lattice-mapping. Via a geometric consideration, we 
have shown conceptually that its inverse exists as long 
as the expansion coefRcients (E^, eA(fc)) of the cluster 
self-energy vanish rapidly in the image-space T^ of the 
projection operator Pij. At closer inspection, the lat- 
tice mapping is thus a two stage process. First, we need 
to determine the expansion coefficients of the cluster self- 
energy. To this end, we will propose a novel interpolation 
technique, which is motivated from the analytical prop- 
erties of the self-energy. The interpolated cluster self- 
energy E^ is then used to compute the inner product 

(Sj: ,ex{kj)) with the eigenfunctions of the projection 
operator Pij, which gives the expansion coefficients of 
the cluster self-energy. Secondly, we need to deconvolute 
the interpolated cluster self-energy on the image space 
If , where we need to determine the optimal value for the 
parameter e. If the latter is too large, the self-consistency 
can not be reached. If e is too small, the lattice-mapping 
will become numerically unstable due to the division of 
small eigenvalues. To solve this problem, we adapt the 
Richardson-Lucy deconvolution algorithm, which inverts 
Eq. ( 23 ) in a numerically stable way. 



s(fc) = ECi?(fc)s^ = E^"''^'VE 



A. Interpolation 

In the context of tight-binding models, one of the most 
successful algorithms to interpolate its band structure is 
the Wannier-interpolation-methocpS. It finds its justifi- 
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c ^ 

K 



pMKy ^ 



The sum over all lattice points can now be split into 
two terms. In the first term, we run over all lattice- 
points within the cluster-radius. In the second term, we 
sum over all the remaining points in the lattice. 



E(fc)=E^"'^"^« 



E e-^^'E^+ E e^^'^'Sfi (28) 

|fl|<-Re \R\>Ra 



If correlations have longer range. Eg. will no longer 
converge on the cluster in real space. This is clearly il- 
lustrated in Fig [t], where we show the self-energy E^ 
for a Nc = 100-site cluster with U/t = 7 for various tem- 
peratures. At high temperatures (T > 0.25), the system 
is only weakly correlated. The self-energy E^ in this 
temperature range is contained within the cluster-radius 
Re = 5. For lower temperatures, it is clear that E^ is 
extends beyond Re- Applying the Wannier interpolation 
scheme according to Eq. (271 to such correlated systems 
is simply not allowed, since the expansion coefficients E^^ 
outside the cluster can not be assumed to be zero. A 



straightforward application of Eq. ( 27 ) will lead to ring- 



ing and eventually to causality violations. The latter was 
observed by Okamata et al.-°', and could only partially be 
resolved by introducing low-pass filtering schemes. The 
applicability of this approach is very limited, due to a 
lack of a general framework to determine these filters. 
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1. Formalism of the interpolation: 

From the previous section, it has become clear that the 
interpolation techniques such as Eq. (27) can only work 
if the function converges on the finite (and often small) 
basis-set. The rate of convergence depends critically on 
the choice of the basis-functions. Consider for example 
the free Green's function G" of the single band Hubbard 
model in Eq. ([3|, 



G°(A?,w) = [ivj + €{k)] 



(29) 



While this Green's function G'^ will converge poorly on 
the cubic-harmonics of the lattice for small frequency -nj, 
it is straightforward to see that [G'^]~^ will be completely 
converged on a 4-site cluster. This simple example shows 
how one can extend the interpolation-idea introduced by 
Okamata et aP^l. Given an injective transformation T, 
we can write 



Hk) = T-'[r[:F]{k) 



K 



(30) 



The method of operation to interpolate a function be- 
comes now clear. Find an injective (and preferably an- 
alytical) transformation T, such that the transformed 
function-values converge on the chosen basis-functions. 
Use this expansion to compute the transformed function- 
values on arbitrary k-points. Finally, apply the inverse 
transformation 7~~^ on the transformed function- values 
in order to obtain the desired interpolated function- 
values on arbitrary k-points. 

This approach has many advantages. First, it pro- 
vides a measure that indicates when the interpolation- 
procedure works or fails. If T[S^] does not converge 
on the basis-set, one is not allowed to perform an in- 
terpolation. Second, this interpolation-procedure does 
not introduce extra information ~ filtering schemes and 
other numerical tricks to ensure causality, on the other 
hand, introduce extra, undesirable structure into the in- 
terpolated functions. By using filtering schemes or other 
numerical tricks to assure causality, we introduce extra 
structure in the function that is to be interpolated, which 
is undesirable. Third, if the transformation T is ana- 
lytical, we will not break the analyticity of the interpo- 
lated function. For Green's functions and their derived 
functions such as the self-energy, analyticity is an impor- 
tant property. In arbitrary interpolation schemes such 
as splines^ or radial-basis expansions'*^, this analyticity 
is often broken. The obtained interpolating function is 
therefore questionable from a physics-point of view. The 
challenge of this approach is naturally the search for a 
correct transformation T- Notice that T can be different 
for different functions, since the only requirements are in- 
jectivity and convergence on the chosen basis-set. In the 



next subsections, we will propose such a transformations 
for the self-energy S. The proposed transformation will 
be motivated by physical and analytical properties of the 
self- energy. 



2. Interpolation on large clusters: 

Since the imaginary part of the self-energy is strictly 
negative in the upper-half of the complex-plane^S 



lm[j:{k,im > 0)] < 0. 



(31) 



we can introduce an injective transformation T that pre- 
serves the analyticity of the self-energySl^ 



r(S) == [S - a i] \ with a > 0. 



(32) 



Due to the property shown in (31 1, the transformation 



T will map the self-energy E into a bounded function, 
irrespective of how spiky the self-energy S is. Notice 
also that we first shift the imaginary part of the self- 
energy down by a i, in order to avoid introducing poles 
due to the Monte Carlo statistical noise. Consequently, 
the function T(S) will now be localized in real space, and 
we can safely perform an expansion of the function 7"(S) 
over cubic harmonics. We have illustrated this process 
in Fig. [8] by applying our interpolation procedure to a 
100-site cluster at a temperature T = 0.2 at half filling. 
In (A) , we show respectively the computed values of the 
cluster self-energy S^ and its interpolation Eg along a 
high-symmetry line in the Brillouin-zone. Notice that 
the imaginary part of the interpolation function remains 
at all times negative! In (B), the transformed function 
7~[E|j] is shown, together with its interpolating function. 
Clearly, the transformation T has reduced the sharp fea- 
tures in the self-energy, and the function has become 
smoother. In (C) and (D), we show the Fourier trans- 
form from respectively the interpolated self-energy T,{k) 
and the transformed values T[Ej^]. The large difference 
in the convergence radii is clear, and shows the effective- 
ness of our indirect approach compared to a direct one. 
This result is not a coincidence. In the appendix, we have 
proven in a rigorous way the point-wise convergence. 



3. Interpolation on small clusters: 

For certain parameter sets, the fermionic sign problem 
prevents the investigation of large enough clusters, for 
which T[E] will converge. In this case, we recommend 
to interpolate the T[E] using cubic splines, instead of in- 
terpolating the latter with the earlier proposed Wannier- 
interpolation. Since T[S] is a much smoother function, 
cubic splines can still perform reasonably well, even in 
the case of small clusters. The self-energy on the other 
hand will not be smooth, and a straightforward spline 
interpolation will lead to overshoots or ringing, which 
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FIG. 8. Interpolation-procedure for the self-energy Ex at the lowest matsubara frequency for a 100-site cluster at a temperature 
T = 0.2, U/t = 7 and t' = at half-filling. (A) The interpolated function E(fc) is a smooth function through the results Ex 
obtained from the QMC cluster solution, where the circles and diamonds represent respectively the real and imaginary part. 
(B) The transformed function T[E] smoothes the self-energy function, making it suitable for a cubic harmonics expansion. (C) 
The Fourier transform of the interpolated function E(fc). Notice that the tails expand much further than the cluster-radius 
Re ~ 5. (D) The Fourier transform of the function T[Ex]. The convergence is reached at Re — 3. 



in turn turn might lead to an acausal self-energy. This 
particular phenomenon has been studied extensively by 
Okamoto et aP^. The ringing might be cured by the use 
of tension splinea^, in which case a tension parameter 
is introduced. It is however important to keep in mind 
that the splines might add extra information into the sys- 
tem, and thus bias the physics. This problem does not 
occur with Wannier interpolation, as long as the Fourier 
coefficients of T[S^] converge on the real space cluster. 



lattice. As a consequence, the lattice-self-energy in the 
DCA breaks the symmetry of the lattice, due to its strict 
parametrization with the coarsegrain patches. The only 
way to resolve this issue in the DCA, is to restrict to 
the few clusters that obey the cluster-symmetry. In or- 
der to remove this undesirable feature in the DCA"*" , we 
symmetrize the self-energy after the interpolation. The 
interpolated cluster-self-energy obeys thus by construc- 
tion the symmetry operations of the lattice. 



^. lattice-symmetry: 



Cluster Deconvolution 



Most of the clusters used in the DCA do generally 
not obey the same symmetry operations as the infinite 



The goal of this section is to present a practical im- 
plementation of the lattice-mapping. As mentioned in 
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the theoretical section of this paper, the lattice map- 
ping is in essen ce t he inversion of the cluster mapping 
defined in Eq. (19). In a common DC A+ calculation, 



we will have much more basis functions than Monte 
Carlo cluster-points. As a consequence, we need to de- 
termine more lattice expansion coefficients than cluster- 
points that are given by the cluster-solver. The inversion 
problem is thus seemingly underdetermined. Therefore, 
we do not attempt to solve Eq. (19) directly, but first 



generalize the coarsegraining equation of the self-energy. 
This is accomplished by rewriting each coarsegraining 
patch as a translation of the patch around the origin, 
i.e 4'ji{k) = 4>^{k — K). Next, we generalize the cluster- 
momentum vector K to an arbitrary momentum vector. 
Using the interpolated cluster self-energy S^ as a substi- 
tute for the cluster self-energy S^ in Eq. (19 1, we obtain 
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Richardson-Lucy iteration 



t{k) 



V 



dk' (j)^{k - k') T.{k') 



(33) 



Any solution of Eq. (33) is thus also a solution of 
Eq. ( 19 ). We should stress that with the exception of the 



continuity of the self-energy, this generalization does not 
introduce any new information as long as the Wannier- 
interpolation converges! With Eq. (33), we have now 



rephrased the lattice-mapping into a deconvolution prob- 
lem. These type of problems are regularly encountered 
in the field of signal theory and image processing and 
various algorithms have been successfully developed to 
address the ill-conditioned deconvolution problenPSl. 



1. The Richardson-Lucy algorithm: 

One of the most common de convo lution algorithms is 
the Richardson-Lucy algorithrnSM! The latter is based 
on a Bayesian inference scheme. Since the patches are 
strictly positive and integrate to unity, we can interprete 
them as a probability distribution function. 



Vfc,fc' : 0g(fc-P) >0, 



1 



Vbz 



BZ 



dk (j)Q{k - k') 



As such, we can apply Bayes theorem and construct a 
conditional probability Q for any given lattice self-energy 



Q{k\k') = 



Mk'-k)j:j{k) 

jg^dk" (j)o(k' -k")Y.(k")' 



(34) 



We should stress at this point that conditional proba- 
bility Q is computed separately for the real and imagi- 
nary part of the self-energy. The conditional probability 
Q{k\K) is then used to construct a new lattice self-energy 
S (fc), given a continuous cluster self-energy S(fc'), 



FIG. 9. Relative error between the cluster self-energy E^ 
and the integrated lattice self-energy f^{K) for the real (open 
symbols) and imaginary (solid symbols) part at 5% doping 
and r = 0.2. 



S(fc) 



dk' Q\k\k')Y.{k'). 



(35) 



BZ 



The idea of the Richardson-Lucy algorithm is now to 
use Eq. (l34| and Eq. (351 in an iterative way. After 



plugging both equations together, we end up with a fixed 
point problem 



S(fc) ^ S(fc) / dk' 



(j)o{k-^k')T.{k') 
J dk" (j)o{k' - k") J:{k" 



(36) 



If the interpolated function S(A:) is now used as our 
initial guess for the lattice self-energy S(fc), Eq. (36) 
provides us with a simple implementation for the 
lattice-mapping. In light of the DCA+ algorithm, the 
Richardson-Lucy deconvolution algorithm has many in- 
teresting properties, that make it an ideal algorithm to be 
used for the deconvolution. First of all, it is a straight- 
forward algorithm that does not need any extra, non- 
physical input . Oth er deconvolution algorithms, such as 
total variatiorP^ESl introduce non-physical penalty fac- 
tors to insure smoothness of the result. Secondly, the 
Richardson-Lucy algorithm conserves the sign of strictly 
positive and negativ e fu nctions. This property can be 
easily proven in Eq. (36), since 4'o{k) is strictly positive. 
Hence, if the initial guess for S(fc) and S(fc') are both pos- 
itive (negative) for all momenta fc, the resulting S(fc) will 
also be positive (negative). Therefore, if the interpolated 
cluster self-energy E(fc) is causal, the lattice self-energy 
will also be a causal function. Third, it has been proven 
that the solution of this iterative scheme converges to 
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FIG. 10. Comparison between the lattice self-energy 

E(fe,7rr), the cluster self-energy T:g{nT) and the coarse- 
grained lattice self-energy at the cluster-momenta T.jf{nT) = 
j:(K,nT) for a 32-site cluster at 5% doping and T = 0.2. 



the maximuni of the likelihood function--. Hence, of 
all lattice self-energies that generate the same cluster 
self-energy after the convolution (coarse-graining), the 
Richardson-Lucy algorithm will produce the lattice self- 
energy that is the most likely to reproduce the cluster 
self-energy. 

Like all other deconvolution algorithms, the 
Richardson-Lucy algorithm is an approximate algo- 
rithm, meaning that the convergence to the exact 
solution is not guaranteed up to an arbitrary precision. 
This is not surprising, since we know that the convolu- 
tion is invertible as long as the expansion coefficients of 
the cluster-self-energy in Eq. (25) decay faster than the 

Consequently, 



eigenvalues of the projection-operator, 
the smaller the cluster, the slower the Richardson-Lucy 
algorithm will converge to a solution and the bigger the 
discrepancy between the coarsegrained lattice self-energy 
f^{K) and the cluster self-energy S^^ obtained from 
the cluster-solver. This phenomenon is illustrated in 
Fig. |9l where we show the relative error in the L2-norm 
between f^{K) and S^. The figure clearly shows that 
the larger cluster converges faster and that the residual 
error between the cluster and coarsegrained self-energy 
decreases with increasing cluster-size. 

In a typical deconvolution, we stop the iteration pro- 
cess if a certain accuracy is obtained. In theory, one could 
use the statistical error of the Monte Carlo integration to 
determine the accuracy. In practice, we have found that 
a relative error below 2.5% delivers in most cases very 
good results. In Fig. 
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we show the lattice self-energy 
for a 32-site cluster, stopped at an accuracy of 2.5%. We 
can clearly observe that the cluster and coarse-grained 
lattice self-energy coincide very well. 



III. APPLICATION: 

A. Convergence of the self-energy and the 
pseudogap: 

One of the most distinctive features of the hole-doped 
cuprates is the emergence of a pseudogap^, i.e. a partial 
suppression of the density of states at the Fermi energy 
at the antinodal points (tt, 0) and (0, tt) in the Brillouin 
zone. This state appears below a temperature T* , which 
rises with decreasing hole doping as the Mott insulating 
half-filled state is approached. The detailed relation be- 
tween the pseudogap and superconductivity remains con- 
troversial. Since superconductivity arises from the pseu- 
dogap state, it is generally believed that understanding 
this unuasual phenomenon is an important prerequisite 
to understanding the pairing mechanism. Recent debate 
has been centered around the question of whether the 
pseudogap is a signature of superconducting fluc tuatio ns 
above T'^SSiSS; ^j. -^^j-^g^^jigj- j^; jg ^ competing phase^^'^. 

Cluster dynamical mean field studies of the single-band 
Hubbard model have found a similar pseudogap opening 
up at the anti nodal po ints at low temperatures in the low 
doping regime^SEfiEMSSl^ j^ these calculations, the pseu- 
dogap originates from a strong momentum-space varia- 
tion of the single-particle self-energy, which, as shown in 
recent DCA calculations by Gull et al.'^Sl, gives rise to 
a momentum-sector-selective metal-insulator transition. 
The DCA+ improves upon the DCA algorithm in that it 
gives a self-energy with smooth and therefore more phys- 
ical momentum dependence, and can therefore provide 
new insight into this problem. In addition, since previ- 
ous studies were limited to relatively small clusters up to 
16 sites, it is important to explore whether the self-energy 
and pseudogap physics is converged on such clusters. 

In Fig. [TT] we plot the imaginary part of the lattice 
self-energy at the smallest Matsubara frequency ujq = ttT 
for various clusters, computed with the DCA (left panel) 
and the DCA"*" (right panel). One immediately observes 
the much more physical smooth momentum dependence 
of the DCA+ results versus the step-function-like nature 
of the DCA results for the self-energy At closer inspec- 
tion, one notices a much more systematic convergence of 
the DCA"*" results with different cluster size and geom- 
etry. While the DCA results for lmE{K) show smaller 
spread at a given ^-point (e.g. at _K^ = (7r,0)), their 
cluster dependence is non-monotonic. In DCA+ , in con- 
trast, |ImE(i4r)| monotonically increases with cluster size 
- a sensible result as longer ranged correlations are sys- 
tematically taken into account. 

Another striking feature of the DCA results is the 
asymmetry for clusters that do not have the full lattice 
symmetry such as the 16B, 20 and 24 site clusters. E.g., 
in the 16B cluster, the asymmetry around (7r/2, 7r/2) as 
one moves along the line from (tt, 0) to (0, tt) is apparent 
and the results in these regions are significantly different 
from those for the symmetric 16A cluster. This asym- 
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FIG. 11. The imaginary part of the lattice self-energy for different clusters at a temperature of T = 0.33 with a hole-doping of 
5% {U/t = 7 and t' /t = —0.15). Two key observations can be made. The DCA"'" produces for all clusters a lattice self-energy 
which follows the lattice symmetry. This is not true in the case of the DCA, which is illustrated in the region of (tt, 0) to 
(0,7r) for the clusters 16-B, 20 and 24. One can also observe that the DCA"'" converges monotonically. The self-energy increases 
systematically with increasing cluster size as longer range correlations are taken into account. This systematic growth of the 
self-energy is harder to detect in the DCA. Therefore, we expect that the DCA^ will lead to a more systematic convergence of 
other physical quantities, such as the pseudo gap transition temperature. 



metry results from the asymmetric arrangenient of the 
two cluster K-points closest to {tt/2, 7r/2) with respect to 
(7r/2,7r/2) (see right hand side of Fig. [I]). This asymme- 
try is completely removed in the DCA~ . 

In addition, with the exception of a small region 
around {tt,tt), the DCA"'" results for the asymmetric 16B 
cluster are almost identical to the results of the fully sym- 
metric 16A cluster. The DCA"'" algorithm restores the 
full lattice symmetry in the results obtained from clus- 
ters that do not have the full symmetry and thus makes 
studies on these clusters much more useful. This, com- 
bined with the improved convergence as a function of 
cluster size allows for much more systematic and precise 
extrapolations to the exact infinite cluster size. 

To further illustrate this point, we now turn to a 
study of the temperature T* below which the pseudo- 
gap starts to form. Here, we define T* as the maxi- 
mum in the temperature dependence of the bulk (q = 
0) magnetic (particle-hole, spin S = 1) susceptibility 
Xphil = 0,T). The downturn in this quantity below 
T* with decreasing temperature signals the suppression 
of low-energy spin excitations, which is also observed in 
experiments to accompany the opening of the pseudo- 
gap in the single-particle spectral weight. In the DCA 
and DCA+ algorithms, Xph is computed from the single 
and two-particle Greens-function G^^ obtained from the 

cluster-solver. Using the notation K = {K,vu), the bare 
two-particle Greens-function Gq{j^ is constructed from a 
pair of interacting cluster Greens functions (for (f = 0) 



while the fully renormalized two-particle Green's func- 
-<ii 

"ph 



tion GU^ is computed as 



gj roi (T1-T2) ^i ^2 (T3-T4) 



Gli{K,K')=U{jyTA 

X Y. (4(^,n)c<,(7?,r2)4,(i?',T3)c<,,(i?',T4)) 



cr.a"'— ± 



The irreducible cluster vertex function Vph{Q = 
0,K,K') is then obtained by inverting the Bethe- 
Salpeter equation on the cluster 



Fp/i — 



'■^Q.ph 



'^ph 



(37) 



where we used a matrix notation in in the cluster mo- 
menta K and K' . The uniform lattice spin susceptibility 
Xph{q — 0) is then calculated from 

Ki,K2 

Here, x" is the coarse-grained bare susceptibility of the 
lattice, 

x\K)= [dkcbK{k)G{k)G{k) 
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FIG. 12. Uniform spin Xph susceptibilities vs temperature 
for different cluster computed in the DCA at 5 percent doping 
{U/t = 7 and t' /t = -0.15). 



FIG. 14. T* versus clustersize computed in the DCA and 
DCA+ at 5 percent doping {U/t = 7 and t'/t = -0.15). 
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FIG. 13. Uniform spin Xvh susceptibilities vs temperature for 
different cluster computed in the DCA^ at 5 percent doping 
(f//i = 7 and i'/i = -0.15). 



the leading correlations are antiferromagnetic and have 
only weak internal i?-dependenceS21, we expect this to be 
a good approximation. 

In Fig. 12 we show results for Xphiq = 0) obtained 
with the DCA for different clusters. One observes a 
strong cluster size dependence and the results are not 
converged even for the largest cluster that can still be 
simulated before the fermonic sign problem begins to 
make the QMC sampling exponentially difficult. The cor- 
responding DCA"*" results are displayed in Fig. 13 Here, 
convergence is reached much sooner. The location of the 
maximum in temperature dependence, T*, is essentially 
independent of the cluster for Nc > 8 (see Fig. 14 1. 
As discussed previously, this directly results from the 
improved convergence of the self-energy in the DCA^ . 
From these results, once the effects of cluster geometry 
are removed in the DCA"*" , it becomes clear that the 
underlying correlations that lead to the pseudogap for- 
mation are short-ranged and well contained in clusters of 



B. Improved fermionic sign-problem 



This procedure to compute the uniform lattice spin 
susceptibility Xph{<l = 0) is the same in the DCA+ as in 
the DCAES. The quantities that enter these equations, 
however, are different between both approaches. In the 
DCA+ , for thermodynamic consistency, one should ap- 
ply the same interpolation procedure to the vertex func- 
tion rpii{K, K') as is done for the self-energy. Here how- 
ever, for the sake of simplicity and in order to focus on the 
effects of the self-energy, we keep the piecewise constant 
dependence oiTph{K, K') that is naturally obtained from 
its extraction from the cluster quantities in Eq. (37) as 
in the DCA. In the S = \ particle-hole channel, where 



The rapidly increasing capability of computers in con- 
junction with the growing sophistication and efficiency 
of quantum Monte Carlo solvers has pushed the lim- 
its of simulations to larger cluster sizes and interaction 
strengths, as well as lower temperatures. As a result, the 
only serious barrier for quantum Monte Carlo calcula- 
tions at low temperatures and away from certain param- 
eter regimes (such as half-filling in the single-band Hub- 
bard model) that remains is the fermionic sign problem^^l^ 
which leads to an exponentially growing statistical error 
with increasing system size and interaction strength, and 
decreasing temperature. 
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The sign problem has posed an insurmountable chal- 
lenge to quantum Monte Carlo calculations of fermionic 
systems, especially for simulations of finite size systems, 
and remains a problem in the DCA approach. The DCA, 
however, was shown to have a less severe sign problem 
than finite size calculations^^, which, in the absence of 
a rigorous mathematical justification, was attributed to 
the action of the mean-field host on the cluster. This has 
enabled simulations of larger clusters at lower temper- 
atures than those accessible with finite size simulations 
and thus has opened new possibilities for gaining insight 
into low temperature phenomena in correlated systems. 

The DCA+ approach is different from the DCA in that 
it generates a more physical self-energy with smooth mo- 
mentum dependence, and the correlations described by 
this self-energy are therefore shorter-ranged than those 
in the DCA. Hence, it is therefore not unreasonable to 
expect a difference in the severity of the sign problem 
between DCA+ and DCA. 

In Fig. 15 we compare the fermionic sign aqmc between 
the DCA and the DCA+ for a 32-site cluster and U = 7t 
for a doping of 5%. At low temperatures, the average sign 
in the DCA+ simulation is significantly larger than that 
of the DCA simulation. As indicated above, we attribute 
this improvement to the smooth momentum dependence 
of the DCA"*" self-energy as compared to the step func- 
tion dependence of the DCA self-energy. From Fourier 
analysis, one knows that the smoothness of a function is 
related to the rate of decay of its Fourier coefficients^. 
More precisely, if a function f is p times differentiable, 
then its Fourier components /„ will decay at least at a 
rate of l/n^"*"^ 



f &CP 



-> 



\fn 



< 






(38) 



Since the DCA+ self-energy has smooth momentum 
dependence and not the step discontinuities of the DCA, 
its Fourier-transform to real space is shorter-ranged than 
that of the DCA and the correlations it describes are 
shorter-ranged. We believe that it is this removal of un- 
physical long-range correlations, which reduces the sign 
problem in the DCA"*" . In any case, with this significant 
reduction in the severity of the sign problem, it is pos- 
sible to study the physics of fermionic systems in even 
larger clusters and at lower temperatures than accessible 
with the DCA. 



IV. SUMMARY AND CONCLUSIONS 

In this paper, we have presented the theoretical 
framework as well as a practical implementation of the 
DCA+ algorithm. It is an extension to the DCA with- 
out the jump discontinuities inherent in the standard 
DCA algorithm that computes a continuous lattice self- 
energy in a self-consistent way. This improvement is 
based on two fundamental differences to the DCA. First, 
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FIG. 15. Temperature dependence of the average fermionic 
sign for Nc = 32 at 5 percent doping {U/t = 7 and t' /t = 
-0.15). 



an explicit distinction is made between the lattice and 
the cluster self-energy. Second, a continuous lattice self- 
energy is determined in a way so that its coarse-grained 
value Ej^ is equal to the cluster self-energy E^ ob- 
tained from the cluster-solver. This constraint makes 
the DCA+ a lgorit hm fundamentally different from previ- 
ous attempta^SHSI to include a continuous self-energy into 
the DCA self-consistency loop that lead to an acausal 
and thus a non-physical self-energy, during the coarse- 
graining of the Greens-function but itself has not been 
coarse-grained. 

The new coarse-graining rules in the cluster-mapping 
of DCA+ require us to reconsider the lattice-mapping in 
the algorithm. As a matter of fact, we have shown that 
a continuous lattice self-energy E(fc) can only be inferred 
from the discrete cluster self-energy E^ if the DCA as- 
sumption of smoothness of the lattice self-energy is satis- 
fied. This has been discussed in the paper using the prop- 
erties of the projection operator Pij that is associated 



with the coarse graining operation in Eq. ( 20 ) . The trans- 



formation of the cluster self-energy into the lattice self- 
energy amounts to inversion of the projection operators 
Pi J. Since this is a singular operator, the lattice map- 
ping is only well-defined as long as the cluster self-energy 
converges on the image-space of the operator, which is 
spanned by the eigenvectors with non-zero eigenvalue. In 
practice the image-space is the space spanned by eigen- 
vectors with an eigenvalue larger than a given parameter 
e. The convergence behavior of the DCA+ algorithm is 
determined by two essential properties of the projectors 
Pij: (1) the dimension of the image-space increases with 
cluster size, which is consistent with the intuitive notion 
that larger cluster can support finer features of the self- 
energy; (2) the delocalization of each eigenvector (r^) 
and the magnitude of its corresponding eigenvalue are 
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anti correlated. Consequently, for large cutoff parameter 
e a more localized cluster self-energy is needed in order 
to have a controlled lattice mapping. Self-consistency in 
the DCA+ can only be reached if the cluster self-energy 
is localized enough to converge on the image space of the 
projection operator. If convergence is not reached, the 
image space of the projector and thus the cluster size will 
have to be increased. Convergence thus provides a use- 
ful measure for the quality of a DCA+ calculation with a 
given cluster. 



Straightforward inversion of the projection operator 
would be numerically imstable, since the projection op- 
erator is a near singular matrix. Thus, in the implemen- 
tation of the lattice mapping in the DCA+ algorithm we 
have followed a different approach, splitting it into two 
numerically stable steps. First, we interpolate the cluster 
self-energy in a controlled way, using an injective trans- 
formation, and next, we deconvolute this interpolated, 
continuous cluster self-energy using the Richardson-Lucy 
algorithm. In both steps convergence within the self- 
consistent loop can be monitored by an objective mea- 
sure. For the interpolation we know that the Fourier 
transform of T[Sj^] — (Sjj — i)^^ has to converge on the 
real-space impurity cluster in order to obtain an accu- 
rate interpolation. For the deconvolution, the difference 
between the coarsegrained lattice self-energy Ej^ and the 
cluster self-energy S^ has to be smaller than the statis- 
tical error of the Monte-Carlo integration. 



To illustrate the benefits of the DCA+ algorithm we 
have investigated the pseudogap phase in a lightly hole- 
doped two-dimensional Hubbard model. Like with the 
DCA, the DCA+ based calculations give a self-energy 
that has strong momentum dependence. However, we 
find that the DCA+ has a much reduced fermionic sign 
problem and thus we can investigate the pseudogap phase 
on larger clusters and in more details than in the standard 
DCA. In the DCA+ the self-energy is continuous in mo- 
mentum space and thus more physical, and it converges 
monotonically and much more systematically with clus- 
ter size than in the DCA. A similarly improved conver- 
gence behavior in the DCA+ is found for the pseudogap 
temperature T* below which the bulk lattice susceptibil- 
ity decreases with decreasing temperature. In the DCA, 
we find that T* has a strong cluster dependence and con- 
verges only for the largest possible cluster sizes. In the 
case of the DCA+ , we observe a much faster convergence 
of T* , which is a direct consequence of the improved con- 
vergence of the self-energy in the DCA+ . From the con- 
vergence property of T* , we can conclude that the corre- 
lations responsible of the pseudogap formation must be 
short ranged and well contained in a cluster size of eight 
sites. This improved convergence in the DCA+ raises the 
hope to do precise extrapolations to the exact infinite 
cluster size limit in future calculations of other proper- 
ties. 
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V. APPENDIX 

Pointwise convergence: Consider a function T and 
an injective, continuous transformation T, such that the 
Fourier components T[J-']r fullfiU, 



Ve>0,3i?eeR: ^ \t[F]r 
\B.\>n.c 



< e 



with T[F]r^ / dke-'^^TiFik)] (39) 



then, 



yk e B, Ve > 0, 3Rc€R: \g{k) - T{k)\ < e 
with .g(fc) =r"i[ ^ exp(ii?fc)r[F]fl, (40) 



R<Rc 



Choose a positive small number e. Since T is a continu- 
ous and invertible function, we know that the T~^ is also 
continuous. Hence, by definition of the this continuity, 
there exists a (5 e IR(j" for this e, such that 



\T[g{k)] - T[Hk)]\ < S ^ \\g{k) - T{k)\ < 



Using the property in Eq. ( 39 ) , we can find a radius 
Re > 0, such that 



E \nFh 

\R\>Ra 



<d. 



(41) 



By the definition of g{k), we have, 

\T[g{k)]-T[Tik)]\ = \ ^ exp(*i?fc)r[^]^ 

-R>|-Rc 

R>Ra 
<S. 



(42) 
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Uniform convergence : Consider a function F 
with the following properties, 

Im[J'(z)] < and \T{z)\ < 6, 6 G R+ 
and an injective, continuous transformation T, 



1 



T{z) = 
such that the Fourier components T[J-']fi fullfill, 



\R\>Rc 



< e 



(43) 



then, 



Ve > 0, 3Rc e R, Vfc e B : \g{k) - J"(fc)| < e 



with g{k) = T^^\ 



R<R^ 



TiFh 



Using condition 
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we have that for all k and for all 



(5 > 0, there exists a critical radius Re such that 

\T[9ik)]-T[Hk)]\ = \ Y, e^^'^T[F]j 
R>\Ro 

R>R^ 
<6. 



Choose the S such that it is smaller then {B + 1)^^/2 
(the minimum convergence length of T~^ on its domain 
is {B + 1)^^). The latter allows us to do an expansion 
around T{k), such that 






E 



9"z 



z=T[T{k)] 



{T[9{k)]-T[Tik)]r 






Hence, we have that for aU k 



\9{k) Hk)\ < E | -^Lfc!i"n+i ('^[g(fe)] - nnkw 

OO 

<E('S+i)"^^'5" 



= iB + ir6j2i{B + l)Sy 



n=0 



iB + l)^S 
l-{B + l)6 



(44) 



The last step is allowed since S < 1/(2 {B + 1)). The 
idea behind the proof is now trivial. Choose an e > 0. 
Pick (5 > such that 



<5 = Min{l/(2(S+l)),-^^^^} (45) 

For this 5, find a radius Re for which, 



E \nF]i 



\B]>Ra 



<5 
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